% This program conducts maximum likelihood estimation of the SVSJ model.
% It is based on the result of "ml_102_3.m" as we use the risk aversion coefficient of that program.
% The jump parameters are from "ml_62_2.m".  Note that we can not use the
% jump parameters estimated in "ml_102_3.m".  This is becase inversion of Y and Z is
% not feasible for some values of jump parameters.
%
% 08-28-05 by Shu Yan.

clear;pack;close all;

% Load the option, return, risk free rate data. For example, opt_1_100 is ATM option with shortest maturity.
% It has the format: opt_1_100=[T r K S C C_b C_a bsV bsV_b bsV_a Vol Opint V X], where C is the
% price, C_b is the bid price, C_a is the ask price, bsV is the BS impv, X is the moneyness ...
load ml_data_5                                    % Results for implied vol errors are in ml_04_YZ_022605

R=ret;
opt_1=opt_1_100(:,[4,3,2,1,5,8]);                 % First option
opt_2=opt_1_105(:,[4,3,2,1,5,8]);                 % Second option
opt_3=opt_2_100(:,[4,3,2,1,5,8]);                 % First option
opt_4=opt_2_105(:,[4,3,2,1,5,8]);                 % Second option
opt_1(:,4)=opt_1(:,4)/365;                        % Change maturity into years
opt_2(:,4)=opt_2(:,4)/365;                        % Change maturity into years
opt_3(:,4)=opt_3(:,4)/365;                        % Change maturity into years
opt_4(:,4)=opt_4(:,4)/365;                        % Change maturity into years
n=length(ret);                                    % length of the sample

% Create a sparse matrix that is used in optimization.  It is used to backed out Y and Z.
Jacb=sparse(1:2:2*n-1,1:2:2*n-1,ones(1,n),2*n,2*n)+sparse(1:2:2*n-1,2:2:2*n,ones(1,n),2*n,2*n)+...
    sparse(2:2:2*n,1:2:2*n-1,ones(1,n),2*n,2*n)+sparse(2:2:2*n,2:2:2*n,ones(1,n),2*n,2*n);

dt=1/52;                                          % Delta t
tau=1/12;                                         % Investment horizon
gam=1.91744;

m_Q=-0.0121556;                              % Estimate from "ml_62_2.m".      
s_Q=0.250408;                                % Estimate from "ml_62_2.m".

m_Q=-0.0154973;s_Q=0.22408;

% Set the options for the likelihood optimization routine.
options=optimset('MaxIter',100,'MaxFunEvals',1000000,...
   'TolX',1e-8,'TolFun',1e-3,'TolCon',1e-6,'DiffMinChange',1e-6,'DiffMaxChange',1);
options=optimset(options,'Largescale','off');
options=optimset(options,'Display','iter');

% Set the options for the pricing error optimization routine.
option_1=optimset('MaxIter',100,'MaxFunEvals',10000,...
   'TolX',1e-6,'TolFun',1e-2,'TolCon',1e-6,'DiffMinChange',1e-6,'DiffMaxChange',1);
option_1=optimset(option_1,'Display','off');
option_1=optimset(option_1,'Display','iter');
option_1=optimset(option_1,'JacobPattern',Jacb);

% parm=[m_Y,k_Y,m_Z,k_Z,s_Y,s_Z,r_12,r_13,r_23] -- Initial guess of the parameters.
parm0=[2.72062  -17.0087  7.02131  -9.59294  0.32398  1.36556  -0.516532  -0.578792  0.181202];

% Bounds for the second group parameters.
LB2=[0.001  0.01 0.2];
UB2=[ 0.12  0.18 1.8];


lb=-10+zeros(1,2*n);ub=10+zeros(1,2*n);

for i=1:20
    disp([i i i i i i i i i i i i i]);
    load ml_102_4_YZ Y0 Z0 parm0
    
    % First use the initial parameters and the options to back out the state variables Y and Z.
    m_Y=parm0(1);k_Y=parm0(2);m_Z=parm0(3);k_Z=parm0(4);
    s_Y=sqrt(parm0(5)^2);s_Z=sqrt(parm0(6)^2);r_12=parm0(7);r_13=parm0(8);r_23=parm0(9);
    
    YZ0=zeros(2*n,1);YZ0(1:2:2*n-1)=Y0;YZ0(2:2:2*n)=Z0;
    [YZ,resnorm,residuals,exitflag,output]=lsqnonlin(@ml_62_4_error,YZ0,lb,ub,option_1,...
        opt_1,opt_2,m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,gam,tau);
    Y=YZ(1:2:2*n-1);Z=YZ(2:2:2*n);
    Y0=Y;Z0=Z;
    parm=zeros(size(parm0));
    
    % Estimate those time series parameters Y and Z.
    [b,bint,u_y,uint,stats]=regress(diff(Y0),[ones(n-1,1),Y0(1:n-1)]*dt);      % We estimate the Y equation by OLS
    parm(1)=b(1);parm(2)=b(2);parm(5)=std(u_y)/sqrt(dt);
    [b,bint,u_z,uint,stats]=regress(diff(Z0),[ones(n-1,1),Z0(1:n-1)]*dt);      % We estimate the Z equation by OLS
    parm(3)=b(1);parm(4)=b(2);parm(6)=std(u_z)/sqrt(dt);
    corr_yz=corrcoef(u_y,u_z);
    parm(9)=corr_yz(1,2);
    
    % Estimate the correlation coefficient between S, Y, and Z.
    parm1=[mean(R)/dt,0.06,0.623];
    [parm2,fval,exitflag,output,Lambda,grad,hessian]=...
        fmincon(@ml_102_3_lkhd,parm1,[],[],[],[],LB2,UB2,[],options,Y0,Z0,R,dt,m_Q,s_Q);
    
    cov_ryz=cov([R(2:n)-parm2(1)*dt,u_y,u_z]);
    parm(7)=cov_ryz(1,2)/parm(5)./(parm2(2)+parm2(3)*mean(Y0))/dt;
    parm(8)=cov_ryz(1,3)/parm(6)./(parm2(2)+parm2(3)*mean(Y0))/dt;
    
    prettymat([repmat(i,size(parm));parm0;parm;parm-parm0]);
    parm0=(parm+parm0)/2;
    save ml_102_4_YZ Y0 Z0 parm0
end

disp('Program is finished!');
disp('-------------------------------------------');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [error,Z_h,parm_h] = ml_62_4_error(YZ,opt_1,opt_2,m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,gam,tau)
% This function calculates the pricing errors in options for SVSJ model.
%
% 04-08-05 by Shu Yan.

n=length(YZ)/2;                                 % Length of Y and Z

Y=YZ(1:2:2*n-1);
Z=YZ(2:2:2*n);

save ml_62_4_YZ_out Y Z

opt=zeros(2*n,6);
opt(1:2:2*n-1,:)=opt_1;
opt(2:2:2*n,:)=opt_2;

Y_s=zeros(2*n,1);Z_s=zeros(2*n,1);
Y_s(1:2:2*n-1)=Y;Y_s(2:2:2*n)=Y;
Z_s(1:2:2*n-1)=Z;Z_s(2:2:2*n)=Z;

T=opt(:,4);
[T_s,I]=sort(T);                                           % Sort so that maturity is in ascending order
[junk,J]=sort(I);                                          % T_s(J)=T

S_s=opt(I,1);K_s=opt(I,2);r_s=opt(I,3);T_s=opt(I,4);C_s=opt(I,5);Impv_s=opt(I,6);
Y_s=Y_s(I);Z_s=Z_s(I);

[pricemodel,Z_h,parm_h]=ml_price(S_s,K_s,r_s,T_s,Y_s,Z_s,m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,gam,tau);
%impvmodel=blsimpv(S_s,K_s,r_s,T_s,pricemodel);
error=C_s-pricemodel;                                      % Pricing error
%error=Impv_s-impvmodel;

error=error(J);
Z_h=Z_h(J);Z_h=Z_h(1:2:2*n-1,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [lkhd,density]=ml_102_3_lkhd(parm,Y,Z,R,dt,m_Q,s_Q)
% This function calculates the likelihood and density for the program.
%
% 08-26-05 by Shu Yan.

% Parameters of the model.
m_S=parm(1);
a=parm(2);b=parm(3);

n=length(R);R_lead=R(2:n);Y_lead=Y(2:n);Y_lag=Y(1:n-1);Z_lead=Z(2:n);Z_lag=Z(1:n-1);
lam=Z_lag.^2;
p=lam*dt;                                      % Jump Probability

% Residuals without jumps.
e_1=R_lead-m_S*dt;

% Covariance matrix of the residuals.
s_11=(a+b*Y_lag).^2*dt;

density=zeros(n-1,1);                     % Initialize density
wt=exp(-p);                               % Initialize weight
% Covariance with the residuals.
XS_iX=e_1.^2./s_11;

density=wt./sqrt(2*pi.*s_11).*exp(-0.5*XS_iX);

num_jps=5;                               % Truncate jump numbers to 5
for j=1:num_jps
    e_1_j=e_1-j*m_Q;
    s_11_j=s_11+j*s_Q.^2;
    
    XS_i_jX=(e_1_j.^2)./s_11_j;
    
    wt=wt.*p/j;                           % Update the weight
    density=density+wt./(2*pi.*s_11_j).*exp(-0.5*XS_i_jX);
end
density=-log(density);

lkhd=sum(density);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%% This program should be first run before the main program %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load ml_data_5                                    % Results for implied vol errors are in ml_04_YZ_022605

R=ret;
opt_1=opt_1_100(:,[4,3,2,1,5,8]);                 % First option
opt_2=opt_1_105(:,[4,3,2,1,5,8]);                 % Second option
opt_3=opt_2_100(:,[4,3,2,1,5,8]);                 % First option
opt_4=opt_2_105(:,[4,3,2,1,5,8]);                 % Second option
opt_1(:,4)=opt_1(:,4)/365;                        % Change maturity into years
opt_2(:,4)=opt_2(:,4)/365;                        % Change maturity into years
opt_3(:,4)=opt_3(:,4)/365;                        % Change maturity into years
opt_4(:,4)=opt_4(:,4)/365;                        % Change maturity into years
n=length(ret);                                    % length of the sample

% Create a sparse matrix that is used in optimization.
Jacb=sparse(1:4*n,ones(1,4*n),ones(1,4*n),4*n,11+2*n)+sparse(1:4*n,1+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4*n,2+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+sparse(1:4*n,3+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4*n,4+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+sparse(1:4*n,5+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4*n,6+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+sparse(1:4*n,7+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4*n,8+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+sparse(1:4*n,9+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4*n,10+ones(1,4*n),ones(1,4*n),4*n,11+2*n)+...
    sparse(1:4:4*n-3,12:2:11+2*n-1,ones(1,n),4*n,11+2*n)+sparse(1:4:4*n-3,13:2:11+2*n,ones(1,n),4*n,11+2*n)+...
    sparse(2:4:4*n-2,12:2:11+2*n-1,ones(1,n),4*n,11+2*n)+sparse(2:4:4*n-2,13:2:11+2*n,ones(1,n),4*n,11+2*n)+...
    sparse(3:4:4*n-1,12:2:11+2*n-1,ones(1,n),4*n,11+2*n)+sparse(3:4:4*n-1,13:2:11+2*n,ones(1,n),4*n,11+2*n)+...
    sparse(4:4:4*n,12:2:11+2*n-1,ones(1,n),4*n,11+2*n)+sparse(4:4:4*n,13:2:11+2*n,ones(1,n),4*n,11+2*n);

dt=1/52;                                          % Delta t
gam=2;                                            % Risk aversion coefficient
tau=1/12;                                         % Investment horizon

% Set the options for the pricing error optimization routine.
option_1=optimset('MaxIter',10000,'MaxFunEvals',1000000,...
   'TolX',1e-6,'TolFun',1e-6,'TolCon',1e-6,'DiffMinChange',1e-6,'DiffMaxChange',1);
option_1=optimset(option_1,'Display','off');
option_1=optimset(option_1,'Display','iter');
option_1=optimset(option_1,'JacobPattern',Jacb);

% parm=[m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,m_S,a,b] -- Initial guess of the parameters.
parm0=[2.97657  -19.2674  6.19024  -4.83268  -0.046278  0.311194  2.71732  0.094411  -0.408095  -0.762775  0.022165];
lb=-100+zeros(1,11+2*n);ub=100+zeros(1,11+2*n);
lb(12:2:11+2*n)=zeros(1,n);

load ml_62_2_YZ_out parm
parm1=parm;

[parm2,resnorm,residuals,exitflag,output]=lsqnonlin(@ml_62_2_error,parm1,lb,ub,option_1,opt_1,opt_2,opt_3,opt_4,gam,tau);

parm=parm2(1:11);
Y=parm2(12:2:10+2*n);Z=parm2(13:2:11+2*n);
    
save ml_62_2_YZ Y Z parm

disp('program is finished!!!');
disp('-------------------------------------------------------------------------------');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%% This program should be first run before the main program %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load ml_data_5                                    % Results for implied vol errors are in ml_04_YZ_022605

R=ret;
opt_1=opt_1_100(:,[4,3,2,1,5,8]);                 % First option
opt_2=opt_1_105(:,[4,3,2,1,5,8]);                 % Second option
opt_3=opt_2_100(:,[4,3,2,1,5,8]);                 % First option
opt_4=opt_2_105(:,[4,3,2,1,5,8]);                 % Second option
opt_1(:,4)=opt_1(:,4)/365;                        % Change maturity into years
opt_2(:,4)=opt_2(:,4)/365;                        % Change maturity into years
opt_3(:,4)=opt_3(:,4)/365;                        % Change maturity into years
opt_4(:,4)=opt_4(:,4)/365;                        % Change maturity into years
n=length(ret);                                    % length of the sample

% Create a sparse matrix that is used in optimization.  It is used to backed out Y and Z.
Jacb=sparse(1:2:2*n-1,1:2:2*n-1,ones(1,n),2*n,2*n)+sparse(1:2:2*n-1,2:2:2*n,ones(1,n),2*n,2*n)+...
    sparse(2:2:2*n,1:2:2*n-1,ones(1,n),2*n,2*n)+sparse(2:2:2*n,2:2:2*n,ones(1,n),2*n,2*n);

dt=1/52;                                          % Delta t
tau=3/12;                                         % Investment horizon
wght=0.1;                                         % Weight of errors in the likelihood function

% Set the options for the likelihood optimization routine.
options=optimset('MaxIter',100,'MaxFunEvals',1000000,...
   'TolX',1e-8,'TolFun',1e-3,'TolCon',1e-6,'DiffMinChange',1e-6,'DiffMaxChange',1);
options=optimset(options,'Largescale','off');
options=optimset(options,'Display','iter');

% Set the options for the pricing error optimization routine.
option_1=optimset('MaxIter',100,'MaxFunEvals',10000,...
   'TolX',1e-6,'TolFun',1e-2,'TolCon',1e-6,'DiffMinChange',1e-6,'DiffMaxChange',1);
option_1=optimset(option_1,'Display','off');
option_1=optimset(option_1,'Display','iter');
option_1=optimset(option_1,'JacobPattern',Jacb);

%m_Q=-0.0121556;                              % Estimate from "ml_62_2.m".      
%s_Q=0.250408;                                % Estimate from "ml_62_2.m".

% parm=[m_Y,k_Y,m_Z,k_Z,s_Y,s_Z,r_12,r_13,r_23,m_Q,s_Q,gam] -- Initial guess of the parameters.
parm0=[2.72062  -17.0087  7.02131  -9.59294  0.32398  1.36556  -0.516532  -0.578792  0.181202 -0.0121556 0.250408 2];

% Bounds for the second group parameters.
LB2=[0.001  0.01 0.2];
UB2=[ 0.12  0.18 1.8];

% Bounds for the jump parameters.
LB3=[-0.1 0.05 1.1];
UB3=[ 0.1  0.5   5];

lb=-100+zeros(1,2*n);ub=100+zeros(1,2*n);

for i=1:10
    disp([i i i i i i i i i i i i i]);
    load ml_102_3_YZ Y0 Z0 parm0
    
    % First use the initial parameters and the options to back out the state variables Y and Z.
    m_Y=parm0(1);k_Y=parm0(2);m_Z=parm0(3);k_Z=parm0(4);
    s_Y=sqrt(parm0(5)^2);s_Z=sqrt(parm0(6)^2);r_12=parm0(7);r_13=parm0(8);r_23=parm0(9);
    m_Q=parm0(10);s_Q=parm0(11);gam=parm0(12);
    
    YZ0=zeros(2*n,1);YZ0(1:2:2*n-1)=Y0;YZ0(2:2:2*n)=Z0;
    [YZ,resnorm,residuals,exitflag,output]=lsqnonlin(@ml_62_4_error,YZ0,lb,ub,option_1,...
        opt_1,opt_2,m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,gam,tau);
    Y=YZ(1:2:2*n-1);Z=YZ(2:2:2*n);
    Y0=Y;Z0=Z;
    parm=zeros(size(parm0));
    
    % Estimate those time series parameters Y and Z.
    [b,bint,u_y,uint,stats]=regress(diff(Y0),[ones(n-1,1),Y0(1:n-1)]*dt);      % We estimate the Y equation by OLS
    parm(1)=b(1);parm(2)=b(2);parm(5)=std(u_y)/sqrt(dt);
    [b,bint,u_z,uint,stats]=regress(diff(Z0),[ones(n-1,1),Z0(1:n-1)]*dt);      % We estimate the Z equation by OLS
    parm(3)=b(1);parm(4)=b(2);parm(6)=std(u_z)/sqrt(dt);
    corr_yz=corrcoef(u_y,u_z);
    parm(9)=corr_yz(1,2);
    
    % Estimate the correlation coefficient between S, Y, and Z.
    parm1=[mean(R)/dt,0.06,0.623];
    [parm2,fval,exitflag,output,Lambda,grad,hessian]=...
        fmincon(@ml_102_3_lkhd,parm1,[],[],[],[],LB2,UB2,[],options,Y0,Z0,R,dt,m_Q,s_Q);
    
    cov_ryz=cov([R(2:n)-parm2(1)*dt,u_y,u_z]);
    parm(7)=cov_ryz(1,2)/parm(5)./(parm2(2)+parm2(3)*mean(Y0))/dt;
    parm(8)=cov_ryz(1,3)/parm(6)./(parm2(2)+parm2(3)*mean(Y0))/dt;
    
    % Estimate the jump parameters and risk aversion coefficient by minimizing pricing errors.
    parm3=[m_Q,s_Q,gam];
    [parm4,resnorm,residuals,exitflag,output]=...
        fmincon(@ml_102_2_error,parm3,[],[],[],[],LB3,UB3,[],options,Y0,Z0,opt_3,opt_4,parm(1:9),tau);
    
    parm(10)=parm4(1);parm(11)=parm4(2);parm(12)=parm4(3);
    prettymat([repmat(i,size(parm));parm0;parm;parm-parm0]);
    parm0=(parm+parm0)/2;
    save ml_102_3_YZ Y0 Z0 parm0
end

disp('Program is finished!');
disp('-------------------------------------------');


